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Abstract 

We discuss two numerical methods, based on a path integral approach 
described in a previous paper (I), for solving the stochastic equations 
underlying the financial markets: the Monte Carlo approach, and the 
Green function deterministic numerical method. Then, we apply the latter 
to some specific financial problems. In particular, we consider the pricing 
of a European option, a zero-coupon bond, a caplet, an American option, 
and a Bermudan swaption. 

1 Introduction 

The evolution law of a financial index, X, is often given by a stochastic differ- 
ential equation, which can be discretized as 

AX = A(X(t),t) At + a(X(t),i) AW, (1) 

where At is the time step, and AW is a Wiener process increment. (In this 
paper, the random variables are denoted by capital letters, and the ordinary 
variables by small ones. Moreover, for the sake of simplicity, we consider only 
the one-dimensional case, but the extension to the multi-dimensional case is 
straightforward.) In a previous paper|l|, hereafter referred as paper I, we have 
seen that the continuous limit of equation ([l]) cannot be written in unambiguous 
form, i.e. it is well defined only if also a discretization rule is given. Therefore, 
from now on, we will always write the underlying stochastic equation in the 
discretized form ([!]), understanding that the continuous limit must be taken. 

Many financial quantities can be defined as the conditional expectation value 
of some functional, <?[JT(t)], of the stochastic process, X(t), obeying the equa- 
tion (|l|). This conditional expectation value can be written as (see, paper I) 



E[g[X(r)} | X(t ) = x ] = 



/+oo p px(t N )=x N r pt N -\ 

dx N Vlaix^y^ir^gixiT^expi- L[x(t), x(t); r]dr \{2) 

-oo J J x(t )=x L Jt ) 

The R.H.S. of Eq. (||) is called path integral[Q, and the functional measure, 
T>[a{x, t) x(r)], means summation on all possible random paths from X(to) = 
xq to X(tisr) — xn, where to, and £jv are the initial and final times, respectively. 
In general, the Lagrangian function, L[x,x;t], is not defined univocally. A 
possible expression is given by 

L[x(t),x(t);t} = — - [x-A(x,t)} 2 . (3) 

Zi (J I JL . I I 

The analytical methods discussed in paper I do not go behind the quadratic 
Lagrangians. This drawback can be partially overcome by using approximate 
analytical techniques such as the perturbative expansion or the saddle point 
approximation, but a more practical and general approach is the numerical one. 

The aim of this paper is to describe two numerical techniques to evaluate 
the path integral above: (a) the Monte Carlo method, which is very general 
and powerful, but has a low precision, and high CPU time requirements; (b) 
the Green function deterministic numerical method (GFDNM), which has been 
recently developed H |L 0, |(|, and it is an advantageous alternative when the 
stochastic process has a low dimensionality. 

In section |[ we recall some general notions on probability theory, and we 
show that the introduction of the conditional expectation value (^) as a path 
integral is just a generalization of the usual concept. In sections ||, and ^, 
we describe the Monte Carlo, and the Green function deterministic numerical 
methods, respectively. Finally, in section ||, some applications of the latter 
method arc discussed. 



2 Preliminary notions on probability theory 

A stochastic process can be defined as a random variable, X (t) , which is function 
of cither a discrete or a continuous time variable, r. The statistical properties of 
the random variable, X(t), at a fixed time, t, are determined by the probability 
density function, or probability distribution function, p(x,t). The stochastic 
process as a whole is characterized by a set of joint probability density functions, 
eventually infinite, 

p(xo,t a ) 
p(xi,ti,x , h) 

p(x 2 ,t 2 ;x 1 ,ti;x ,to) (4) 

with ti € [to,ijv]j which must satisfy the Kolmogorov compatibility conditions, 
J p(x N ,t N ; . . .;xi,ti;. . . ;x ,to) dxi = 

p(xN,tN] • ■ ■ ; aj»+i,ti+i;xi_i,ti_i; . . . ; xq, h) . (5) 
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The conditional probability density functions are defined as the ratios 

p(x N ,t N ;...;x i ,t i | Xi-i,ti-i;...',x ,t ) = — . (6) 

p(Xi-!, ti-i; . ..;xo,to) 

Let us now assume that the time variable is discrete. We can define the 
expectation value of a function of the stochastic process, g(X(t^), . . . , X(to)), 
as 

E[g(X(t N ),...,X(t ))] = 

dx N ...dx g(x N , ...,x ) p(x N ,t N ; ...;x , t Q ) , (7) 



and the conditional expectation value as 
E[g(X(t N ), . . . , X(t )) | X(ti-i) = Xi-i; . . . ; X(t ) = x ] = 

dx N ...dxi g(x N , ...,x ) p(x N ,t N ; ...;Xi,ti \ a^i,^^; . . .;x ,t ). (8) 



The equation above defines a function of the ordinary variables, Xi-\, . . . ,xq. If 
these variables are substituted by the random variables, X(ti-i), . . . , X (to), the 
conditional expectation value (||) becomes a function of the stochastic process. 
In this case we will use the shorter notation 

E\g(X(t N ), . . . , X(to)) I x(u-i); x(t ) } = 

E[g{X{t N ), X(to)) | X(ti-i) = X(ti-i)i ■ ■ -;X(to) = X(to)]. (9) 

Finally, we can define the expectation value with fixed initial and final points 

(x N , t N | g(X(t N ), . . . , X(t )) | x , t ) = 

dxN-i ...dxi g(xN, ■ ■ ■ ,x ) p{xN,t N ] ...;xi,tx \ x , to). (10) 

Some relations among the quantities above are 
E\g(X(t N ),...,X(t ))] = 

dxi-i ...dx E[g(X(t N ), . . . , X(t )) \ Xfc-i) = Xi_i; . . . ; X(t ) = x ] 

x p(x i _i,t i _ 1 ; . . .;x ,t ) , (11) 

and 

E\ g (X(t N ),...,X(t ))\X(t ) = x ] = 

dx N {x N ,t N | g(x N , ...,x ) I x , t ) ■ (12) 



Let us now consider two particular types of stochastic processes: the Markov 
process, and the Gaussian process. 
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The Markov process is defined as a stochastic process such that 

p(xi,ti | ari_i,ti_i; . . .;x ,t ) = p(xi,ti | (13) 

i.e. the conditional probability density function at the time, U, depends only 
on the next earlier time, and not on the whole previous history of the process. 
It then follows that the conditional probability density with fixed initial point 
can be written as 

N 

p(x N , t N ;...;xi,t! | x , t ) = JJ p{x u U \ (14) 

The function, p(xi,ti | Xi-i,ti-i), is also called transition probability. The 
following properties of a Markov process can be easily proved 

E[E[ gi (X(t N ),...,X(U)) \X(t i )}g 2 (X(t i ),...,X(t ))\X(t ) = x ] = 

E[ 9l (X(t N ), . . . , X{U)) g 2 (X(U), . . . , X(t )) | X(t ) = x ], (15) 

and 

E[g(X(t N ), . . . , X(U)) | X(U) = Xi ; X(t) = x] = 

E[g(X(t N ),...,X(t i ))\X(t i )=x i ], (16) 

for every t < U. 

Finally, a stochastic process is called Gaussian process if all its joint proba- 
bility density functions are multivariate Gaussian distributions. 

2.1 Expectation value of a functional of random paths 

If the time variable, r, is continuous, instead of a function of a stochastic process 
we should more properly speak of a functional, g[X(r)] (note the square brackets 
in the notation), of a random path. Its expectation value and its conditional 
expectation values can be written as a generalization of the equations in the 
previous section, i.e. 

E [g[ x ( T )}} = lim /••■/ T\dx l g(x N ,...,x ) p(x N ,t N ;...;x ,t ), (17) 

J J i=0 



(x N ,t N I g[x(r)} | x ,t ) 

N — l 

= lim /. . ./ TT dxi g(x N , . . . ,x ) p{x N ,t N ; . . . ;xi,ti \ x Q ,t ) (18) 

N^oo J J JLX 

= ^ V[x(r)] g[x(r)] p[x(r)l (19) 

J J x(t )=x a 
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and 



/+oo 
dx N (x N ,t N | g[x(r)\ | x Q ,t ) 
-OC 



(20) 



Moreover, for a Markov process, we have a straightforward generalization of 
Eqs. (H), and @, i.e. 



E 



E 



gi[X(r) | r G Mat]] | *(i<) g 2 [A(r) |rG [to, *i] ] I *(*(>)= ar 



5i[X(T)|rG[t i) i JV ]] ff2 [X(T)|Te[to ) ti]]|X(to)=a;o , (21) 



and 



S [X(t) | r G [ti,t N ] } | X(ti) = x, ; X(t) = x 



E 



g[X(j) | re [ti,t w ]] | X(ti)=Xi], 



(22) 



where t < f,-, and, for the sake of clarity, we have explicitly written the intervals 
of variation of the time variable, r. 

The expression (|l^) represents a definition of the path integral but 
the latter is only a formal way to write the limit of the former discretized 
expression. However, for a Wiener process, and more generally, for a Langevin 
process, described in the next section, an exact mathematical measure can be 
defined (Wiener integral) 0. 

The functions g(xjj, . . . , xq), and p{x^,tj^; . . . ;xi,ti | xo,to) represent a 
discretization of the functionals, 5[x(r)], and p[x(r)]. We want to stress here 
that, in general, both the explicit form of the functionals and the discretization 
procedure are not defined univocally. A more detailed discussion of this problem 
can be found in paper I. 



2.2 Wiener and Langevin processes 

A stochastic process can be Gaussian and Markovian at the same time. The 
Wiener process, for example, is defined by the initial probability density func- 
tion, p(x, to) = S(x), and the Gaussian transition probability, 

9 {x lM | Xt-uU-i) = - 7 == ex p{- ( TaS } } ' (23) 

where At = U — t,_i. 

Another Markov, in general non-Gaussian, process is the solution of the 
Langevin equation (|l|). An explicit expression for the joint probability den- 
sity function corresponding to this equation does not always exist. However, a 
general expression of the transition probability for small time steps (short-time 
transition probability), correct up to O(Ai), is given by (see, paper I) 

p(xi,U | Xi-i,U-i) ~ 

1 / (xj - Xj-i - ife-i.tj-i) At) 2 

y/2ira{x t - 1 ,t t ^ 1 ) 2 At ^ \ 2 oix^xM-i) 2 At 
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Therefore, 



p(x N ,t N ; ...;xi,ti | x ,t ) ~ (2Tra(x i ^ 1 ,t i ^ 1 ) 2 At) N/2 

..... -Xi-i -^ \- l)Af)2 l . (25) 
' 7 J 2 f j( a ; l _i,t i _i) 2 At [ v y 




Note that the expression (|24|) is equal to (|23|) when A = 0, and er = const. 
In this case the short-time transition probability ( p4| ) is exact, and the Wiener 
process is the solution of the Langevin equation. 

From now on, if not explicitly specified, we will consider only Langevin 
processes. The probability density functional for a Langevin process, apart 
from a normalization factor usually included into the measure, can be written 
as t 

p[x{r)} ~ exp j- jf N L[x{t), x(t); r]dr| ; (26) 

then the Eq. ( p0| ) becomes equal to the conditional expectation value (||). We 
recall that the Lagrangian (||) , with the discretization rule (|2^) , corresponds to 
the pre-point formulation of the path integral (see, paper I). 



3 The Monte Carlo method 

Monte Carlo is a technique for the numerical computation of mathematical 
quantities using random numbers ||. Its basic idea lies on two important limit 
theorems of probability theory: the law of large numbers, and the central limit 
theorem. As a result of these theorems, if f{X) is a function of a random 
variable, X, and the random numbers, x^ r \ are sampled from its probability 
distribution, for large M the average 



M 

= I7 E^ (r) )' ( 27 ) 



is an estimator of /i = E[f(X)]. Moreover, even if the random variable has a 
quite general probability distribution, the averages obtained by different sam- 
plings are distributed according to a normal probability density function with 



M 



mean /x, and variance a = E 
given by the quantity, 



Finally, an estimation of o is 



- — ^r~n^ iJir ' } " l)2 - m 



„2 

" M(M - 1) 



These theorems can be generalized to the case of multivariate probability 
distribution functions. 
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3.1 Monte Carlo integration 

The aim of this introductory section is to describe the use of the Monte Carlo 
method to perform the numerical integration of an ordinary function. We briefly 
discuss both the plain and the importance sampling methods. 

3.1.1 Plain integration 

Let us consider the ordinary integral, 

rb 



J a 



x) dx. (29) 



If a^ 1 ), . . . ,a;( M ) are M random numbers uniformly distributed in the interval 
[a, b], as a result of the limit theorems, we have 



M 

I = E u [f(X)]^I=-J2f(x {r) ) (30) 

r=l 

(the indices, u, means that the expectation value is computed on the uniform 
probability distribution) , and the statistical error is of the order of the standard 
deviation, 

M 

(31) 



3.1.2 Importance sampling 

We must often calculate an integral of the form, 

I=( f{x)p{x)dx, (32) 

J a 

where p(x) is some probability density function defined in the interval [a, b). It 
can happen that in the most part of the interval, the function, p(x), is almost 
zero. In this case, with the uniform sampling defined above, a large number of 
points gives a negligible contribution. In order to overcome this limit, we can 
use the importance sampling method: if x^\ . . . , x^ are random numbers 
sampled from the distribution function, p(x), for the limit theorems again, we 
have 



1 M 

I = E p [f(X)}^I=-J2^ x{r) ^ ( 33 ) 



r=l 

and the statistical error, 



A/ 
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3.2 Random number generation 

The previous section has shown that, in general, we need random numbers 
sampled from an arbitrary probability distribution. Uniform random number 
generators are implemented on almost all computers, but the sampling from 
more complex probability distributions must be performed by appropriate algo- 
rithms. In this section we describe two fundamental algorithms of this kind. 

3.2.1 Acceptance-rejection algorithm 

This method was proposed by von Neumann^] to obtain random numbers with 
a given probability distribution, p(x). Let p(x) be defined in a finite interval 
[a, b], and bounded by a value, U. Then random numbers with probability 
distribution, p(x), can be obtained by the following algorithm 

i. Draw x from a uniform distribution in [a, b]. 

ii. Draw p from a uniform distribution in [0, U]. 

iii. If p < p(x) then accept x else reject it, and goto i. 

In order to have a better efficiency, the upper bound, U, must be as close as 
possible to the maximum of p{x). Unfortunately many points must be rejected, 
and the method is not very efficient, especially in the case of multivariate prob- 
ability distributions. Moreover, it can be applied only to bounded distributions 
with a finite range. 

3.2.2 Metropolis algorithm 

A more efficient and general method is given by the Metropolis algorithmpO|. 
This is related to an important property of the Markov processes: after a number 
of time steps large enough, the final probability density of a stochastic system, 
which evolves according to a discrete time Markov process, is independent on 
the number of steps, and on the initial configuration. The Metropolis algorithm 
is a solution of the inverse problem: is it possible to construct a Markov process 
which, after a number of steps large enough, yields a configuration with a given 
probability density function, p(x)7 The solution is not unique, but the Metropo- 
lis procedure is particularly simple to implement. It consists of the following 
steps: 

i. Start at the initial time from a point x. 

ii. Draw a random number, q, from the uniform distribution in the inter- 
val [0, 1], and compute x' — x + D(2q — 1), where D is some arbitrary 
parameter. 

iii. If p(x') > p(x) then goto v. 

iv. Draw another random number, p, from the uniform distribution in the 
interval [0, 1]. If p(x') < p p{x) then put x' = x. 
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v. Put x' in the list of random numbers, rename x' to x, and goto ii. 

The random numbers, a;', have the required probability density function. Actu- 
ally, a certain number of steps must be performed until convergence is attained, 
i.e. a certain number of initial random points must be rejected. The parameter, 
D, is arbitrary, but it has been empirically proved that an appropriate value 
of D should give a ratio of acceptances between 50% and 70%. If we have a 
multivariate distribution, and hence a multi-dimensional array of points, each 
point should be moved in turn. The new array corresponds to one step of the 
Markov process. 



3.3 Expectation value computation 

A continuous stochastic process, X(t), can be specified either by its evolution 
law (stochastic differential equation), or by its probability density functional 
(path integral formulation). Thus, it can be approximated either by the solution 
of the discretized Langevin equation (|l|), or by the discretized stochastic process 
with the multivariate probability density function (p5|). 

According to the limit theorems, we can write the following approximate 
expression for the conditional expectation value ( |20| ) , 

M , (r) (t-)n 

E[g[X(r)] | X(t ) = x ] - km £ ^ ' " ■ *° ^ (35) 

M— >oc * — * 1V1 

r=l 

where the random paths, . . . , x%\ are selected either by solving Eq. (|l|), 
or by knding the stochastic process with the probability density (E5I); in both 

(r) 

cases with a kxed starting point, x y Q ' — x . We will now briefly describe how 
these two problems can be solved by the Monte Carlo methods. 



3.3.1 Langevin equation approach 

A solution, Xq , . . . of the Langevin equation (Q) can be obtained by the 

following algorithm: 

(r) 

i. Put i = 0, and Xq = xq. 

ii. Draw a random number, z, from the normal distribution with mean and 
variance At. 

iii. Take ^ = xf ] + A{xf\u) At + a(x<f \ t t ) z. 

(r) 

iv. Store xW x . 

v. If i < (N — 1) then put i = i + 1, and goto ii. 



An approximation of the conditional expectation value (20) is given by the 
summation in Eq. (|35|), where the sum extends to all paths generated. 



9 



3.3.2 Path integral approach 

The discretized conditional expectation value, defined by Eqs. (|l8|), and (|2^), is 
simply a multi-dimensional integral. Therefore, we can use the Metropolis algo- 
rithm to find an ensemble of discretized paths with the multivariate probability 
density function ( p5f ) |Tl| . An approximation of the conditional expectation value 
( p0| ) is given again by the sum in Eq. ( J35| ) (see also Ref. ]20[). 

We want to point out here that the Markov process underlying the Metro- 
polis algorithm must not be confused with the Markov process describing the 
physical process, and obeying the Langevin equation (Q). The former is only 
a formal device to obtain the required joint probability density. A sample 
path corresponding to the first stochastic process is given by the set of vec- 
tors, (xq 1 ^, . . . , ) , . . . , {xq \ . . . , Xjf'), while a sample path for the second 

(r) (r) 

process is simply given by the single array, Xq , . . . , x N . 



4 Green function deterministic numerical me- 
thod 

The main drawback of the Monte Carlo method is the CPU time requirement. 
We must remind that to gain one order of magnitude in the precision we need to 
increase the CPU time of two orders. In this section we describe an alternative 
method |j| which has some advantages in low dimensional problems. In the 
following we will call indifferently Green function or transition probability the 
conditional probability density (|l3|), since actually this function represents the 
Green function of the partial derivative equation corresponding to the Langevin 
equation (Q). 



4.1 Expectation value computation 

Let us consider a general Markov process. The probability density functional, 
p[a;(T)], can be approximated by the product of transition probabilities given in 
Eq. (p"4|). If we assume that the functional, g[x(r)], can be discretized in the 
form, 

N 

g[x(T)]~l[gW( Xi ), (36) 

i=0 

which includes most of the interesting cases, then an approximation of the con- 
ditional expectation value (EOh is given by 



N N 



tj-l), 



E[ 9 [X{t)] I X(t ) = x ] * f ... ( f[d Xl g (N) (x N ) f[p{xj,tj | x^, 

J %=i j=i 

(37) 

where the function, p, is 

p(xi,U | Xi-uti-i) = p{xi,U | Xi-x,ti-i) g^Hxi-x). (38) 
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Let us now consider a single integration 

J dxi p(x i+1 ,t i+1 | Xi,U) p(xi,ti | Xi-i,U-i) . (39) 

If we approximate this integral by using a numerical quadrature rule, we obtain 
the following algebraic relation 

M 

£ a# ~p^ 1] *h > (4°) 

7=1 

where the matrices, p W, are defined by 

PaP = P(Za,U+l | £/3,*i), (41) 

the quantities, w a , and z a , are the weights and the grid points, respectively, 
associated with the integration rule, and a, /?, 7 = 1, . . . , AT. In conclusion, the 
expression ([37]) can be written as 

M 

E[g[X(r)} I X(t ) = z a ] c £ gW cg-^ . . . , 

(42) 

where G^l = w a p^, and 5 7 ^' ) = g( N \z 7N ). Therefore, we have reduced 
the evaluation of the expectation value of a functional to the product of N 
matrices with dimension M. By starting the calculation from the left, we need 
to memorize just linear arrays, while the matrix elements, G a p, can be computed 
step by step. In practice, we adopt the following algorithm: 

i. Put u a = gi N] (a = 1, . . . , M), and i = N - 1. 

M 

ii. Put v a = up Gf a (a = 1, . . . , M). 

0=1 

iii. If i > then put u a — v a (a = 1, . . . , A'/), i = i — 1, and goto ii. 
Here the arrays, u Q , and u OJ are two working vectors. 



4.1.1 Path dependent options 



In sections 5.4 
tions 



and 5.5 we will discuss two examples of path dependent op- 



The pricing of this type of derivative securities requires to evaluate the 
expectation value of functionals containing some constraints. In general, their 
explicit expression is very involved, and an exact analytical treatment is not 
possible. Instead a pure numerical approach implies only a small difference in 
the algorithm described above. In particular, for the cases considered in this 
paper, it becomes simply 
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i. Put u a = g a ' {a = 1, . . . , M), and i = N — 1. 

I M 

ii. Put Wa = f(z a ,ti), and v a — max I it^ G 




) 



(a = l,...,M). 

iii. If z > then put u a — v a (a — 1, . . . , M), i = i — 1, and goto ii. 

The only difference is in the point ii, where we have now a test operation, and 
the function, f(z a ,ti), depends on the problem considered. 

4.2 Computational details 

The method described above is quite general. Let us now consider the particular 
case of a stochastic process obeying the Langevin equation ([j]) . Here we can use 
the approximate expression (|24|) for the short-time transition probability. 

Since the interval of integration in Eq. ( ^ ) is infinite, the numerical inte- 
gration should be performed by a quadrature rule for improper integrals (for 
instance, the Gaussian quadrature). On the other hand, since the integrand 
is essentially a narrow Gaussian of width a \f~Ki whose central position, xi, 
moves on the whole interval, we need a grid dense enough to give an accurate 
quadrature everywhere. Therefore we are forced to take a uniform distribu- 
tion of grid points with a lattice spacing, Ax — Xi — x^\ ~ er V At; where 
ef = min t)), and the interval, X, is the finite range of integration 

xEX, T6[t ,tjv] 

due to the finite number of points. As a result the stochastic process is confined 
in a box, but an interval large enough gives negligible corrections. A good choice 
for the quadrature formula is the trapezoidal rule, which yields very accurate 
results with Gaussian functions, and, in general, with functions which are zero, 
with all their derivatives, out of some range. 

Actually, the relation between At and Ax, and the need of a finite range of 
integration are two connected problems: if we fix the number of grid points and 
the interval, I, then the value of At is fixed by the relation At ~ Ax 2 /a 2 (in 
practice, a ratio W 2 At/ Ax 2 = 1 gives an accuracy greater than 1%, and already 
with 1.25 we get at least 8 digits). In other words, if we take At going to zero, 
the Gaussian part of the transition probability becomes strongly peaked, and 
we need a very large number of points to obtain a good precision. This means 
that we cannot take At as small as possible, and the systematic error which 
depends on At can be significant (this systematic error must not be confused 
with the numerical error in the quadrature, which depends only on the ratio 



As a final remark, we note that the definition above for a can be sometimes 
source of troubles. In particular, when the function, a(x,t), has some zeros or 
a very wide range of variation, the time step, At, can be very large. In this case 
a larger number of points than usual is necessary. 



a 2 At /Ax 2 ). 
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For financial problems, we typically have 50-100 time steps, and about 50 
grid points. Therefore, a one-dimensional problem is reduced to the computation 
of the product of 100 matrices of dimension 50, which can be handled on a PC 
in a few seconds. The result is usually accurate to a level of 10~ 3 without any 
particular caution. If a very high precision is required we can either improve 
the approximation of the short-time transition probability or decrease the time 
step, At. 

In the first case, we must expand the short-time transition probability in 
powers of At and Ax, by requiring that this expression satisfies up to any given 
order the partial derivative (Fokker-Planck) equation associated to the Langevin 
equation (|l|) [0, [| . The expansion coefficients are analytical but rather cumber- 
some. Each coefficient increases the precision of about one order of magnitude, 
but there is a trade-off between the required precision and the complexity of the 
analytical expressions. 

In the second case, we must either increase the number of grid points ac- 
cording to the relation between At and Ax, or expand the short-time transition 
probability over some basis of interpolating functions^). In the latter method 
the relation a 2 At ~ Ax 2 is not essential anymore. Obviously we have an ad- 
ditional cost, in terms of CPU time requirements, for the computation of the 
expansion coefficients. Again there is a trade-off between the required precision 
and the complexity of the method. 

4.3 Advantages and limits 

The main advantages of the GFDNM with respect to the Monte Carlo methods 
are: 

- higher velocity and accuracy; 

- the path dependent derivative securities are easily handled; 

- the solution for all initial values of the stochastic variable is obtained in a 
single iteration. 

The drawback of this approach appears in the multi-dimensional case. In a 
d-dimensional problem, in general, we need matrices of dimension M d . For 
example, in four dimensions with 50 grid points, we have matrices of dimension 
50 4 = 6250000. Although such matrices are sparse (as a result of the Gaussian 
form of the transition probability), i.e. a large number of their elements are 
zero, they must be stored on the hard disk, increasing the computing time. 

However, we do not need the transition probability matrix itself, but just 
to multiply this matrix with a vector of M d points, which, at least in four 
dimension, can be easily stored in the RAM of a PC. On the other hand, the 
short-time Green function is an analytical function which can be computed 
wasting some CPU time, but without troubles of memory requirement. 

Obviously, as d increases, the dimension of the vector M d grows up, and the 
problem becomes quickly intractable. In this case, the only numerical method 
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available is the Monte Carlo, which always gives an answer (but not necessarily 
correct). The comment in bracket is not sarcastic: in high dimensional problems 
the random sampling allows to get a definite answer, but the configurations 
taken into account are only a very small part of the possible ones. Thus, the 
probability of missing some critical set of paths is not negligible. 



5 Numerical results 

In the following we apply the Green function deterministic numerical method 
to some specific problems. 

First, in order to compare the numerical results to the exact ones, we consider 
three examples with analytical solutions. In particular, we compute the price 
of a European option on a non-dividend-paying stock in the Black and Scholcs 
model] 21 1, and the prices of a zero-coupon bond, and a caplet in the Vasicek 
model j§2|j. 

Second, we apply the method to two path dependent derivatives. In partic- 
ular, we compute the prices of an American option on a non-dividend-paying 
stock in the Black and Scholes model, and of a Bermudan swaption in the Va- 
sicek model. 

We want to stress that the choice of the above financial models is only due to 
the fact that we can compare the numerical results to the exact solutions. The 
numerical code, of course, is very general, and does not depend on this choice. 
The program also assume a possible time dependence of the parameters, and 
could be optimized for constant parameters saving much CPU time. 



5.1 European option 

A European option gives the holder the right (and not the obligation) to buy 
or to sell an underlying asset on a certain date (exercise date, or maturity) 
for a certain price (exercise price). In the first case it is called call option, in 
the second one put option. If S(T) is the price of the underlying asset at the 
maturity, T, and x is the exercise price, the expected value of a European put 
option at the time, t, in a risk-neutral world is 

Oe {a u t, T) = E \e~ r (T -^ max( X - S(T), 0) | S(t) = s t ] , (43) 

where the risk-free rate of interest, r, is assumed constant for the whole life of 
the option. The functional, g[s(r)], is then given by 

g[s(r)} = e~ r ^ max( X - s(T),0). (44) 

The expression above can be discretized simply as 

N-l 

g[a(r)] * max( X - s N ,0) J[ e~ rAt , (45) 
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and the linear array, g a , is 

ff f)=max(x-s a) 0). (46) 

Furthermore, since the coefficients of the Lagrangian for the Black and Sc- 
holes model are <j(S) — a S, and A(S) = fiS, the matrix, G a p, is given by 

n 1 / (s a - Sf} -iJ,Sf) At) 2 \ 

G af 3 = w a = exp i r At S , (47) 

where s Q are the grid points, w a are the weights of the integration rule, and 
a, ft = 1, . . . , M. In Table |l] we show a comparison between analytical and 
numerical results. 



5.2 Zero-coupon bond 



A zero-coupon bond is a contract which yields a certain amount [principal) on a 
certain date {maturity) in the future. If, for the sake of simplicity, the principal 
is equal 1, the price at the time, t, is given by the conditional expectation value 
of the functional (see, paper I) 



rT 

r i m ~ I r ( T ) dT 
g[r(r)J = e J * 



(48) 



where T is the maturity. The functional above can be discretized by the pre- 
point rule as 



N-l 



g[r(r)} ~ J] 



-j-j At 



(49) 



i=0 



with At 



T-t 

AT 



Hence the linear array, g& , is 



/ 1 \ 
1 

7 W _ 



(50) 



V 1 / 



The coefficients of the Lagrangian (||) are cr(r) = a, and A(r) = a(6 — 7-). 
Therefore the matrix, G a p, is given by 



G a = w c 



1 



V2vro- 2 At 



exp 



(fa - - a(b - zp) At) 2 
2 a 2 At 



zpAt\, (51) 



where the grid points. Now, by performing the product (ft2j), we ob- 

tain directly the zero-coupon bond price, P(z a , t, T), for each initial short term 
interest rate, z a . 

In Table | we show a comparison between analytical and numerical results. 
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5.3 Caplet 



A more complex example is the computation of the caplet price. An interest 
rate cap is a contract which guarantees that the rate charged on a loan does not 
exceed a specified value, the cap rate. The cap can be viewed as a portfolio of 
European put options on zero-coupon bonds. The individual options comprising 
a cap are referred to as caplets. The expected value of a caplet is given by the 
expectation value of the functional (see also paper I) , 

g[r(r)} = e~ £ r{r)dT ( X - P(r T , T, a)) o(x - P(r T , T, «)) (52) 

= e -J t r ^ dT max (x- P(r T ,T,s),0). (53) 

In order to calculate the value of the caplet, we first need the price of the zero- 
coupon bond at the time, T, which can be computed by the method described 
above. This procedure (unlike a Monte Carlo approach) gives directly the price 
for any grid point, z a , and allows direct integration over this variable. Once the 
zero-coupon bond price has been computed, the vector, g&\ is given by 

5 W=max(x-P(^,T,s),0), (54) 

while the matrix, G a p, is still given by the expression (jsjj). 

In Fig. Ill we show a comparison between analytical and numerical results. 



5.4 American option 

A path dependent option is an option whose value depends on the past history 
of the underlying asset, not just on its value on exercise. As a first example of 
a path dependent option, we consider an American put option, i.e. an options 
which can be exercised at any time up to the expiration date. Unfortunately 
an exact analytical expression for the price of an American option does not 
exist. Therefore the importance of a powerful numerical technique which is 
able to handle this problem is evident. The Monte Carlo method is not very 
appropriate in this case, since has the disadvantage that early exercise features 
are difficult to implement. The techniques usually adopted are based on the 
binomial trees |23| or on the resolution of partial differential equations by finite 
differences p4| . The Green function deterministic numerical method described in 
the previous sections is particularly efficient, in terms of convergence properties, 
memory requirements, accuracy, and implementation. 



In section 5.1 we have seen that, if S(T) is the price of the underlying asset 
at the maturity, T, and x is t ne exercise price, the expected value of a European 
put option at the time, t, in a risk-neutral world is given by Eq. (|43]). On the 
other hand, if the option is American, its value cannot be defined so easily. 
From a computational point of view, however, it can be obtained as follows. 
The option price, in general, is evaluated by starting at the final time, T, where 
the values of the European and the American options coincide, and are simply 
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given by 

O a (S(T), T, T) = E (S(T),T, T) = max( x - S(T), 0), (55) 

and working backward in time steps, At. If the option is American, we also 
need to check at any time, and for any value of the stock price, whether early 
exercise is preferable to holding the option for a further time step. Therefore 
the procedure for the calculation is the following: 

i. We compute the value of the option at the time, T — At, if it is not 
exercised: 

E [e- rAt A (S(T), T, T) \ S(T - At) = s T _ At } (56) 

(we recall that the risk-free rate of interest, r, is constant here). 

ii. We compute the value of the option at the time, T — At, if it is exercised: 

max(x- s T _ A t,0). (57) 

iii. The correct expected value of the option at the time, T — At, is given by 

Oa (st-auT- At,T) = 

max(£ [e- rAt A (S(T),T,T) | S(T - At) = s T _ Af ] , 

max(x - ST—At ,0)j, (58) 

i.e. we check if it is convenient to exercise or not the option. 

iv. We put T = T — At, and we iterate the procedure until the initial time. 

In conclusion, by recalling the algorithm and the notations given in section pi. l.l| , 
we have 

fl W = max( X -s Q ,0), (59) 
= max( X -s Q ,0), (60) 

and the matrix, G a p, is the same of Eq. ([f7j). 

In Table | we show a comparison of the numerical results obtained by a 
finite-difference method j24|, a binomial method ]23[|, and the Green function 
deterministic numerical method. 

5.5 Bermudan swaption 

A swap contract is an agreement between two companies to exchange cash flows 
in the future according to a defined formula. The simplest example of swap is 
the plain vanilla interest rate swap. In this case a first party agrees to pay to 
a second party cash flows equal to interest at a predeterminate fixed rate on a 



17 



principal for a number of years. At the same time the second party agrees to 
pay to the first party cash flows equal to interest at a floating rate on the same 
principal for the same period of time, in the same currency. Actually, there 
is an enormous number of different swap types that can be invented. For the 
sake of simplicity, we will consider only the plain vanilla. However, the method 
described in this section can be easily extended to more complex cases. 

A swaption is an option on a swap contract. It gives the holder the right to 
enter into or to terminate a swap contract at a certain time in the future. In 
the following we will consider a Bcrmudan swaption, which is a particular non- 
standard American option. In a Bcrmudan swaption early exercise is restricted 
to certain dates during the life of the swap, usually the reset dates. 

Let us then consider a plain vanilla interest rate swap settled in arrears, 
with a principal, Q, a fixed rate, K, the floating rate equal to the London 
Interbank Offer Rate (LIBOR), the fixed rate payment dates, if-, . . . ,tf/, and 
the floating rate payment dates, tf,... , tjy . We can assume that the floating 
base rate underlying the swap is the appropriate rate to use for discounting. 
This is a common assumption (see, for example, section 5.2 in Ref . |p3| ) , and 
considerably simplifies the valuation procedure. Therefore the floating rate, I/j, 
at the time, tf, is given by 



L t ^L(R(tf),tf,tf +1 ) = -L 



1 



P(R(tf),tf,tf +1 ) 



1 



(61) 



The swap price, W(t), at the time, t, is given by the sum of the actualized 
values of all remaining payments after t. If we define the integers, Nk, Nl, 
with 1 < N K < N K , and 1 < N L < N L , such that fc£ , and tk- are the dates 
of the first fixed rate and floating rate payments after t, respectively, we have 



N L N K 

w{t)= ]T wf(t,tf)- wftttf), 



(62) 



i=N L 



where 



W % K (t,tf) = E 



E 



,-/_' R(r)dr 
i-t K 

- _* RMdr 
» Jt 



L<_i A L | R(t) 



Vr < t 



Q K A K | R(t) =r T , Vr < t 



, (63) 



(64) 



A K =tf - tti, Al = tf - tf_ v and t$ = tf - A K , iff = - A L . Note 



that the floating rate payment at the time tf is made according the rate at the 
beginning of the period, i.e. in arrears. Since R(t) is a Markov process, the 
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expressions above are equivalent to 
E 



E 



- h' R ^ dT QLi-iA L | R(t) = r T ; R(U-i)= r t< _ 1 



(65) 



Wf[i,t?)=E 



- R r) di 

3 J* V ' 



(66) 



By exploiting the expression ( |6l| ) for the LIBOR, and the properties (pl|), and 
@, the Eqs. ©, and ©, become 

Wf (*,*f) = 



- / J fl(r)dr 



E 



P{R{tt x ),ttirf) 



I R® = rj 



t<t 



i-l 



fl(r)dr | R{1)=r _ 



E 



(/T 



P(i) 



if_,<i<*f 



(i'(r F ,t,tf_ 1 )-P(r f ,t,tf)) 



-P(rj,ltf) 



(67) 



(68) 



In conclusion, the calculation of the expected value of a swap contract can be 
reduced to that of zero-coupon bonds, and can be made by using the procedure 



described in section 5.2 



Let us now consider a Bermudan swaption which gives the right to terminate 
the swap contract defined above at various dates, Ti, . . . , TV. If we assume that 
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the exercise dates, T, coincide with some floating rate payment dates, tf, as 



usual, the swap price at the time, t = T, depends on the short term interest 
rate, rj, only, i.e. W(t) = W(rj,t). The expected value at the time, t < Ti, of 
the corresponding European swaption is simply 



(r t ,t,,T N )=E 



- «(r) dr ^ {wmTN)iTN) } 0) | R{t) = n 



(69) 

In the case of a Bermudan swaption, we also need to check at all times, Tj, and 
for any possible value of the short term interest rate, whether early exercise is 
preferable to holding the swaption for a further time interval. The procedure 
is analogous to that of an American option. The main difference is that the 
underlying asset price is not the stochastic variable itself, but the swap price. 
Therefore the swap price must be calculated in advance for all exercise dates of 



the swaption. We use an algorithm similar to that of section 4.1.1 with 

gW = max(VK(r Q ,T w ), 0), (70) 
u;W = max(W(r a ,TJ,Q), (71) 

and the matrix, G a p, given in Eq. (pl|). The only difference is that the test in 
step ii is made only at the exercise dates of the swaption. 

In Table ^, we show the numerical results obtained for a plain vanilla interest 
rate swap in the Vasicek model compared with the analytical ones. In Table § 
and Table ^ we show a comparison of the numerical solutions obtained for the 
European and the Bermudan swaptions, respectively, with the results obtained 
by a semi-analytical computation. 

5.6 Greeks 

Any financial institution has the problem of hedging the risk of its portfolio 
of derivative securities. For this reason it needs to know the sensitivity of 
the portfolio to the changes of the underlying asset prices, the time, and the 
market conditions. This sensitivity is usually measured by calculating five hedge 
parameters (greeks): 

• delta, is the rate of change of the derivative security price with respect to 
the price of the underlying asset; 

• gamma, is the rate of change of the portfolio's delta with respect to the 
price of the underlying asset; 

• theta, is the rate of change of the value of the portfolio with respect to 
time; 

• vega, is the rate of change of the value of the portfolio with respect to the 
volatility of the underlying asset; 

• rho, is the rate of change of the value of the portfolio with respect to the 
interest rate. 
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The Green function deterministic method gives directly the price of a derivative 
security as a discretized function of the initial value of the underlying variable, 
and the whole discrete time evolution of the price is calculated step by step. 
Moreover, the price variation with respect to other parameters can be obtained 
by changing these parameters, and performing a new computation. Therefore, 
the calculation of the greeks does not add any further complication. 

5.7 Conclusions 

We have described two numerical methods, based on the path integral formu- 
lation given in Ref. p|], to calculate the conditional expectation value of a 
general functional: the Monte Carlo method, and the Green function determin- 
istic numerical method. Moreover, we have shown some practical applications 
of the latter to the pricing of derivative securities. In order to compare the 
analytical and the numerical results, we used two solvable financial models: the 
Black-Scholes, and the Vasicek models. However, the GFDNM works also in 
more complex cases, when an analytical solution does not exist. The numerical 
results are very accurate, and can be even improved by using some particular 
techniques §, % §. The method has been applied to one-dimensional problems, 
but the extension to the d-dimensional (d ~ 2 — 4) ones is straightforward. An- 
other important feature of the GFDNM is that it gives directly the conditional 
expectation value of a functional for all initial values of the stochastic variable. 
Finally, the case of the path dependent derivative securities can be handled with 
only small changes in the code. 

In conclusion, the GFDNM is a very powerful technique for low dimensional 
problems, as usually the derivative security pricing. On the other hand, the 
Monte Carlo method is the only possible, although very slow and imprecise, for 
high dimensional problems. 
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Table 1: Comparison of the analytical Black-Scholes solution with the numerical 
solutions obtained by a binomial method, and the Green function deterministic 
numerical method for a European put option with x = 10, t = yrs, and 
T = 0.5 yrs. 



Stock price Analytical Binomial[23| GFDNM 



6.0 3.558 3.557 3.557 
8.0 1.918 1.917 1.917 
10.0 0.870 0.866 0.871 
12.0 0.348 0.351 0.349 
14.0 0.128 0.128 0.129 

Note. The Black-Scholes parameters are r = 0.1, and <r = 0.4. The number of grid points 
for the GFDNM is 201, and the CPU time is about 3 seconds on a Pentium 133. 



Table 2: Comparison of the analytical Vasicek solution with the numerical so- 
lution obtained by the Green function deterministic numerical method for a 
zero-coupon bond with t = yrs, and T — 0.5 yrs. 



Short term 
interest rate 


Analytical 


GFDNM 


Relative error 


0.02 


0.9884 


0.9886 


2 x 10~ 4 


0.04 


0.9797 


0.9797 


8 x 10~ 5 


0.06 


0.9710 


0.9709 


8 x 10~ 5 


0.08 


0.9625 


0.9622 


3 x 10" 4 


Note. The Vasicek parameters are a = 
points is 51, and the CPU time is much lc 


0.5, b = 0.05, and cr 
ss than 1 second on a 


= 0.03. The number of grid 
Pentium 133. 



Table 3: Comparison of the numerical Black-Scholes solutions obtained by a 
finite-difference method, a binomial method, and the Green function determin- 
istic numerical method for an American put option with \ — 10, t = yrs, and 
T = 0.5 yrs. 



Stock price Finite difference §J Binomial j§3| GFDNM 
(TO 4000 1000 4.000 
8.0 2.095 2.096 2.093 
10.0 0.921 0.920 0.922 
12.0 0.362 0.365 0.364 
14.0 0.132 0.133 0.133 



Note. The Black-Scholes parameters are r = 0.1, and a = 0.4. The number of grid points 
for the GFDNM is 201, and the CPU time is about 3 seconds on a Pentium 133. 
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Table 4: Comparison of the analytical Vasicek solution with the numerical solu- 
tion obtained by the Green function deterministic numerical method for a swap 



with Q = 1, K = 0.045, t = yrs, = t% = 10 yrs, Ajf = Al = 0.5 yrs, and 
N K = N L = 10. 



Short term 
interest rate 


Analytical 


GFDNM 


Relative error 


0.02 


0.01064 


0.01066 


2 x 10~ 3 


0.04 


0.01037 


0.01037 


4 x 10~ 4 


0.06 


0.01010 


0.01009 


9 x 10~ 4 


0.08 


0.00984 


0.00982 


2 x 10" 3 



Note. The Vasicek parameters are a = 0.5, b = 0.05, and a = 0.03. The number of grid 
points is 51, and the CPU time is about 3 seconds on a Pentium 133. 



Table 5: Comparison of a semi-analytical Vasicek solution with the numerical 
solution obtained by the Green function deterministic numerical method for a 
European swaption with Q = 1, K = 0.045, t — yrs, t$ = t$ = 10 yrs, 
Ak = Al = 0.5 yrs, and Nk = Nl = 10, which gives the right to terminate 
the swap just after the 8 th payment date. 





Semi-analytical 


GFDNM (51 pts) 


GFDNM (101 pts) 


Short term 


Swaption 


Swaption 


Relative 


Swaption 


Relative 


interest rate 


price 


price 


error 


price 


error 


0.02 


0.00409 


0.00423 


3 x 10" 2 


0.00412 


7 x 10" 3 


0.04 


0.00393 


0.00405 


3 x 10~ 2 


0.00396 


8 x 10~ 3 


0.06 


0.00377 


0.00389 


3 x 10~ 2 


0.00380 


8 x 10~ 3 


0.08 


0.00362 


0.00372 


3 x 10" 2 


0.00365 


8 x 10" 3 



Note. The Vasicek parameters are a = 0.5, b = 0.05, and <r = 0.03. The CPU time with 51 
grid points is about 3 seconds on a Pentium 133. 



Table 6: Comparison of a semi-analytical Vasicek solution with the numerical 
solution obtained by the Green function deterministic numerical method for a 
Bcrmudan swaption with Q = 1, K = 0.045, t = yrs, t$ = t$ = 10 yrs, 
Ak = A^ = 0.5 yrs, and = = 10, which gives the right to terminate 
the swap at the time iff, and just after the 2 nd , 4 th , 6 th , and 8 th payment dates. 





Semi-analytical 


GFDNM (51 pts) 


GFDNM (101 pts) 


Short term 


Swaption 


Swaption 


Relative 


Swaption 


Relative 


interest rate 


price 


price 


error 


price 


error 


0.02 


0.01559 


0.01609 


3 x 10" 2 


0.01569 


6 x 10" 3 


0.04 


0.01494 


0.01539 


3 x 10~ 2 


0.01503 


6 x 10~ 3 


0.06 


0.01431 


0.01472 


3 x 10~ 2 


0.01440 


6 x 10~ 3 


0.08 


0.01371 


0.01408 


3 x 10~ 2 


0.01380 


7 x 10" 3 



Note. The Vasicek parameters are a = 0.5, b = 0.05, and <r = 0.03. The CPU time with 51 
grid points is about 3 seconds on a Pentium 133. 
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Figure 1: Comparison of the analytical Vasicek solution (solid line) with the 
numerical solution (dots) obtained by the Green function deterministic numer- 
ical method for a caplet with \ = 1-001, t = yrs, T — 0.5 yrs, and s = 1 yrs. 
The Vasicek parameters arc a = 0.5, b = 0.05, and a = 0.03. The number of 
grid points is 51, the relative error is of the order of 10~ 2 , and the CPU time is 
much less than 1 second on a Pentium 133. 
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